## R code for analyses
rm(list=ls())
setwd("C:/Users/your_work_space")
getwd()

## Packages:
#***********

library(vegan)
library(adespatial) 
library(ggplot2)
library(forecast)   
library(spacodiR)
library(ade4)
library(rethinking) ## for Bayesian credibility interval
library(mediation)
library(dplyr)

############################
## Phenological trap data ##
############################

#----------
## FLOWERS:
#----------
FF0 = as.data.frame(read.table("SuppInfo_Appendix_S6.txt",h=T))
head(FF0) ; dim(FF0)

## Filter species with Nr of basket >= N in at least one year:
# (NOT SURE I NEED TO DO THIS)
# - - - - - - - - - - - - - - - - - - - - - - - - - - - - - -
# The filter that is applied is consistent with the methodology of Wright & Calderon (2006) Ecology Letters
# for BCI (also a 50 ha plot with 200 originally placed phenology baskets [50 were added in gaps later]).
# We subset the community to species that were found in at least 10 phenology collection baskets in at least
# one year (irrespective of phenophase - i.e. either fruits of flowers). Chang Yang (2016) Jounal of Ecology
# applied a similar filter, but with a more conservative number of basket (5 or more), because the plot of Fushan
# is smaller. Applying the filter yielded 129 species for flower records, 367 for fruiting records, and 391 in total. 

A2 = FF0 %>% filter(QUANTITY >= 0) %>% group_by(CODIGO,YR,TRAP) %>% tally(wt = NULL) ; A2
A3 = as.data.frame(A2) ; head(A3) ; dim(A3)
head(A3) ; dim(A3)

nrsp   = length(table(A3$CODIGO)[table(A3$CODIGO)>0]) ; nrsp
namesp = names (table(A3$CODIGO)[table(A3$CODIGO)>0]) ; namesp

N = 10 # Min Nr of traps in one year or more
i=10;j=1
specs = c()
for(i in 1:nrsp){
  sub  = A3[A3$CODIGO%in%namesp[i],] ; sub
  nyrs = length(table(sub$YR)[table(sub$YR)>0]) ; nyrs
  yrs1 = names (table(sub$YR)[table(sub$YR)>0]) ; yrs1
  max  = 0
  for(j in 1:nyrs){
    sub2 = sub[sub$YR%in%yrs1[j],] ; sub2
    if(nrow(sub2)>=N){max=c(max,1)}
  }
  if(sum(max)>=1){specs = c(specs,namesp[i])}
}
specs 

## Now create new matrix of flower data with only the species
## found in at least 10 traps in one year or more:
FF <- FF0[FF0$CODIGO%in%specs,] 
head(FF); dim(FF) ; dim(FF0)
names(table(FF$CODIGO)[table(FF$CODIGO)>0]) ## verification (= length(specs))
names(table(FF$YR)[table(FF$YR)>0]) ## verification

## Check if Min Number of traps >= N:
verif=c() 
for(u in 1:length(specs)){
  samp2    = FF[FF$CODIGO%in%specs[u],]
  verif[u] = length(table(samp2$TRAP)[table(samp2$TRAP)>0])
}
min(verif) ## Must be >= 10

#--------------------------------------------
## Same procedure for the fruit observations:
#--------------------------------------------
TT0 = as.data.frame(read.table("SuppInfo_Appendix_S7.txt",h=T))
TT0$spf = NA ## add empty column for "seed per fruit" 
TT0$sdm = NA ## add empty column for "seed dry mass"
head(TT0) ; dim(TT0)

## Filter species with Nr of basket >= N in at least one year:
# (NOT SURE I NEED TO DO THIS)
# - - - - - - - - - - - - - - - - - - - - - - - - - - - - - -
A2 = TT0 %>% filter(QUANTITY >= 0) %>% group_by(CODIGO,YR,TRAP) %>% tally(wt = NULL) ; A2
A3 = as.data.frame(A2) ; head(A3) ; dim(A3)
head(A3) ; dim(A3)

nrsp   = length(table(A3$CODIGO)[table(A3$CODIGO)>0]) ; nrsp
namesp = names (table(A3$CODIGO)[table(A3$CODIGO)>0]) ; namesp

N = 10 # Min Nr of traps in one year or more
i=10;j=1
specs = c()
for(i in 1:nrsp){
  sub  = A3[A3$CODIGO%in%namesp[i],] ; sub
  nyrs = length(table(sub$YR)[table(sub$YR)>0]) ; nyrs
  yrs1 = names (table(sub$YR)[table(sub$YR)>0]) ; yrs1
  max  = 0
  for(j in 1:nyrs){
    sub2 = sub[sub$YR%in%yrs1[j],] ; sub2
    if(nrow(sub2)>=N){max=c(max,1)}
  }
  if(sum(max)>=1){specs = c(specs,namesp[i])}
}
specs 

TT <- TT0[TT0$CODIGO%in%specs,] 
head(TT); dim(TT) ; dim(TT0)
names(table(TT$CODIGO)[table(TT$CODIGO)>0]) ## verification (=length(specs))
names(table(TT$YR)[table(TT$YR)>0]) ## verification

verif=c()
for(u in 1:length(specs)){
  samp2 = TT[TT$CODIGO%in%specs[u],]
  verif[u] = length(table(samp2$TRAP)[table(samp2$TRAP)>0])
}
min(verif) ## Must be >= 10


## Appendix S1: List of species, their family and life form guild:
# - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - -
head(TT); dim(TT)
sp1 = names(table(TT$CODIGO)[table(TT$CODIGO)>0]) ; sp1
nsp = length(sp1) ; nsp

apps1 = as.data.frame(matrix(nrow=nsp,ncol=3))
colnames(apps1) = c("species","family","guild")
apps1$species = sp1

i=1
for(i in 1:nrow(apps1)){
  wf              = which(TT$CODIGO%in%apps1$species[i])[1]
  apps1$family[i] = TT$fam[wf]
  apps1$guild[i]  = TT$guild[wf]
}

head(apps1);dim(apps1)

#write.table(apps1,"SuppInfo_Appendix_S1.xls")

#-------------------------------------------------------
## Seed Dry Mass (SDM) and Seed(s) Per Fruit (SPF) data:
#-------------------------------------------------------

SPF.SDM = as.data.frame(read.table("SuppInfo_Appendix_S2.txt",h=T))
head(SPF.SDM) ; dim(SPF.SDM)

## Nr and names of species in the trait matrix
nrsp = length(row.names(SPF.SDM)) ; nrsp
nams = row.names(SPF.SDM) ; nams

## Fill spf and sdm colums in TT, using correspondence with the SPF.SDM file:
i=1
for(i in 1:nrsp){
  w1 = which(TT$CODIGO%in%nams[i],)
  if(length(w1)>0){
    TT$spf[w1] = SPF.SDM$SPF[i]
    TT$sdm[w1] = SPF.SDM$SDM[i]
  }
}
head(TT); dim(TT)
TT[1:100,]
length(which(is.na(TT$sdm)))/nrow(TT) ## no missing values remaining
length(which(is.na(TT$spf)))/nrow(TT) ## no missing values remaining

#-----------------------------------------------------------------------
## Calculate the Nr of plots in which each flowering species is observed
## during each census
#-----------------------------------------------------------------------

head(FF);dim(FF)
Nplots = length(table(FF$YRMON)[table(FF$YRMON)>0]) ; Nplots
Plots  = names(table(FF$YRMON)[table(FF$YRMON)>0]) ; Plots

FP    = as.data.frame(matrix(ncol=4,nrow=Nplots))
colnames(FP) = c("YR","MON","YRMON","FP")
FP$YRMON = as.numeric(Plots)
FP$YR    = floor(as.numeric(Plots))
FP$MON   = as.numeric(Plots)-floor(as.numeric(Plots))
FP$MON   = FP$MON/0.08
head(FP) ; dim(FP)

barplot(12-table(FP$YR),ylab="Nr Missing Months")
w=which(table(FP$YR)==12);w;length(w)
years=as.numeric(names(table(FP$YR))[w]);years ## years with complete flower data 

FP=FP[which(FP$YR%in%years),]
head(FP);dim(FP)
Plots2 = names(table(FP$YRMON)[table(FP$YRMON)>0]) ; Plots2
Nplots2 = length(Plots2) ; Nplots2

nrsp = length(names(table(FF$CODIGO)[table(FF$CODIGO)>0])) ; nrsp
nams = names(table(FF$CODIGO)[table(FF$CODIGO)>0]) ; nams

## Matrix that will contain flower production data 
SP.FP = as.data.frame(matrix(0,ncol=nrsp,nrow=length(Plots2)))
colnames(SP.FP) = nams
rownames(SP.FP) = Plots2
head(SP.FP[,1:10]) ; dim(SP.FP)

for(i in 1:nrow(FP)){ ## For each census...
  
  subF = FF[FF$YRMON==Plots2[i],] # dim(subF)

  nams2 = names(table(subF$CODIGO)[table(subF$CODIGO)>0])
  nrsp2 = length(nams2)
  
  for(j in 1:nrsp2){ ## For each species in the census i...
  
  subFF  = subF[subF$CODIGO%in%nams2[j],] ; dim(subFF)
  ntraps = nrow(subFF)
  
  w99 = which(colnames(SP.FP)%in%nams2[j])
  if(length(w99)>0){ SP.FP[i,w99] = ntraps }
  
  }
}

head(SP.FP) ; dim(SP.FP)
length(which(colSums(SP.FP)>0))
SP.FP = SP.FP[,which(colSums(SP.FP,na.rm=TRUE)>0)]
SP.FP = SP.FP[,rev(order(colSums(SP.FP,na.rm=TRUE)))]
plot(colSums(SP.FP,na.rm=TRUE))

head(SP.FP[,1:10]) ; dim(SP.FP)
min(rowSums(SP.FP,na.rm=TRUE))
min(colSums(SP.FP,na.rm=TRUE))

#----------------------------------------------------------------------------
## Create a table with total seed mass for each species for each plot:
#----------------------------------------------------------------------------
head(TT);dim(TT)  ## fruit data with only the species present in >= 10 traps
                  ## in one year or more (generated just above)
                  ## and with sdm and spf data

Nplots = length(table(TT$YRMON)[table(TT$YRMON)>0]) ; Nplots
Plots  = names(table(TT$YRMON)[table(TT$YRMON)>0]) ; Plots

SRM    = as.data.frame(matrix(ncol=3,nrow=Nplots))
colnames(SRM) = c("YR","MON","YRMON")

SRM$YRMON = as.numeric(Plots)
SRM$YR    = floor(as.numeric(Plots))
SRM$MON   = as.numeric(Plots)-floor(as.numeric(Plots))
SRM$MON   = SRM$MON/0.08
head(SRM) ; dim(SRM)

nrsp = length(names(table(TT$CODIGO)[table(TT$CODIGO)>0])) ; nrsp
nams = names(table(TT$CODIGO)[table(TT$CODIGO)>0]) ; nams

table(SRM$YR)
w=which(table(SRM$YR)==12);w;length(w)
years=as.numeric(names(table(SRM$YR))[w]);years

SRM9 = SRM[SRM$YR%in%years,] ## New SRM with only complete years for fruit data
head(SRM9);dim(SRM9)

Plots2  = names(table(SRM9$YRMON)[table(SRM9$YRMON)>0]) ; Plots2
Nplots2 = length(Plots2) ; Nplots2

SP.SRM = as.data.frame(matrix(0,ncol=nrsp,nrow=length(Plots2)))
colnames(SP.SRM) = nams
rownames(SP.SRM) = Plots2
head(SP.SRM[,1:10]) ; dim(SP.SRM)

for(i in 1:Nplots2){ ## For each YRMON census...
  
  subTT = TT[TT$YRMON%in%Plots2[i],] 
  nams2 = names(table(subTT$CODIGO)[table(subTT$CODIGO)>0])
  nrsp2 = length(nams2)
  
  for(j in 1:nrsp2){ ## For each species in the census i...
    
    subTTT = subTT[subTT$CODIGO%in%nams2[j],] 
    w1  = which(subTTT$PART==1)  ; Q1  = 0
    w2  = which(subTTT$PART==2)  ; Q2  = 0
    w3  = which(subTTT$PART==3)  ; Q3  = 0
    w7  = which(subTTT$PART==7)  ; Q7  = 0
    w10 = which(subTTT$PART==10) ; Q10 = 0
    w11 = which(subTTT$PART==11) ; Q11 = 0
    
    if(length(w1)>0) { subT1  = subTTT[w1,]  }
    if(length(w2)>0) { subT2  = subTTT[w2,]  }
    if(length(w3)>0) { subT3  = subTTT[w3,]  }
    if(length(w7)>0) { subT7  = subTTT[w7,]  }
    if(length(w10)>0){ subT10 = subTTT[w10,] }
    if(length(w11)>0){ subT11 = subTTT[w11,] }
    
    if(length(w1)>0) { Q1  = sum(subT1$QUANTITY *subT1$spf *subT1$sdm)  }
    if(length(w2)>0) { Q2  = sum(subT2$QUANTITY            *subT2$sdm)  }
    if(length(w3)>0) { Q3  = sum(subT3$QUANTITY *subT3$spf *subT3$sdm)  }
    if(length(w7)>0) { Q7  = sum(subT7$QUANTITY            *subT7$sdm)  }
    if(length(w10)>0){ Q10 = sum(subT10$QUANTITY*subT10$spf*subT10$sdm) }
    if(length(w11)>0){ Q11 = sum(subT11$QUANTITY           *subT11$sdm) }
    
    w99 = which(colnames(SP.SRM)%in%nams2[j])
    if(length(w99)>0){ SP.SRM[i,w99] = sum(c(Q1,Q2,Q3,Q7,Q10,Q11)) }
    
  }
  
}

head(SRM9);dim(SRM9) 
head(SP.SRM) ; dim(SP.SRM)
length(which(colSums(SP.SRM)>0))
SP.SRM = SP.SRM[,which(colSums(SP.SRM,na.rm=TRUE)>0)]
SP.SRM = SP.SRM[,rev(order(colSums(SP.SRM,na.rm=TRUE)))]

head(SP.SRM[,1:10]) ; dim(SP.SRM)
min(rowSums(SP.SRM,na.rm=TRUE))
min(colSums(SP.SRM,na.rm=TRUE))

#--------------------------------------------------------------------
## Calculate community-level flower production (FP2) for each census:
#--------------------------------------------------------------------

log1 = log(SP.FP+1) ; log1                        ## log of the flower production data created above
SP.FPA = decostand(log1,"standardize",na.rm=TRUE) ## z-score transform data for each species
head(SP.FPA[,1:10]) ; dim(SP.FPA)

## Standardized log-transformed community-level flower production:
SP.FPB = SP.FPA 
FP2 = as.numeric(rowSums(SP.FPB,na.rm=TRUE)) ; FP2 

#--------------------------------------------------------------------------------
## Idem for fruit (same way & for flowers) => Calculate community-level 
## seed production (SRM2):
#--------------------------------------------------------------------------------

log2 = log(SP.SRM+1)
SP.SRMA = decostand(log2,"standardize",na.rm=TRUE) 
head(SP.SRMA[,1:10]) ; dim(SP.SRMA)

## Standardised log-transformed community-level seed production:
SP.SRMB = SP.SRMA 
SRM2 = as.numeric(rowSums(SP.SRMB,na.rm=TRUE)) ; SRM2

#****************************************************************
## ANALYZE TRENDS in total seed production at the community level
#****************************************************************

# - - - - - - - - - - - - - - - - - - - - - - - - - - -
## STL and inter-annual trends (not used at the moment):
# - - - - - - - - - - - - - - - - - - - - - - - - - - -
my.vector3 <- ts(1:length(SRM2), frequency = 12, start = c(2000, 1)) 
my.vector3[1:length(SRM2)] <- SRM2
my.vector3

stl1 <- stl(my.vector3, "per") ; plot(stl1)
season <- as.numeric(stl1$time.series[,1]) ; plot(season,type="l")
intera <- as.numeric(stl1$time.series[,2]) ; plot(intera,type="l")

# - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - -
## B-splines and linear regression with quadratic approximation:
# - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - -
D = as.data.frame(cbind(SRM9$YRMO,SRM2)) 
colnames(D) = c("X","Y")
head(D);dim(D)

sd.Y    <- sd(D$Y,na.rm=TRUE)
means.Y <- mean(D$Y,na.rm=TRUE)
D$Y     <- as.numeric(decostand(D$Y,"standardize",na.rm=TRUE))
#forecast::InvBoxCox((D$Y*sd.Y+means.Y),lambda=lamb, biasadj = FALSE, fvar = NULL)
plot(stats::density(D$Y))
hist(D$Y) 

## (1) B-splines:
# - - - - - - - -
num_knots <- 20
knots_list <- quantile(D$X,probs=seq(0,1,length.out = num_knots))
Q <- bs(D$X,knots=knots_list[-c(1,num_knots)],degree=2,intercept = TRUE)
Y <- as.numeric(D$Y) ; mean(Y,na.rm=TRUE) ; sd(Y,na.rm=TRUE)
  
m4.7 <- quap(
  alist(
  Y ~ dnorm(mu,sigma), mu <- a + Q %*% w, a ~ dnorm(6,10), w ~ dnorm(0,1), sigma ~ dexp(1)
), data=list(D=D,B=Q),
start=list(w=rep(0,ncol(Q))))
  
summary(m4.7)
post <- extract.samples(m4.7)
w <- apply( post$w , 2 , mean )
mu <- link( m4.7 )
mu_PI <- apply(mu,2,PI,0.95)
precis(m4.7)
  
mu1     <- mu_PI
#mu1[1,] <- InvBoxCox((mu_PI[1,]*sd.Y+means.Y),lamb, biasadj = FALSE, fvar = NULL)
#mu1[2,] <- InvBoxCox((mu_PI[2,]*sd.Y+means.Y),lamb, biasadj = FALSE, fvar = NULL)
  
plot(0:12,col="white",
     ylim=c(min(D$Y,na.rm=TRUE)*1.2,max(D$Y,na.rm=TRUE)*1.2), 
     xlim=c(2000,2022),col.axis="white",xlab="Year",ylab="Total Seed Mass")
axis(side=1, at=c(2000:2022), col="black",labels=c(2000:2022)) 
axis(side=2, at=c(-4,-3,-2,-1,0,1,2,3,4), 
     col="black",labels=c(-4,-3,-2,-1,0,1,2,3,4))

points(D$Y~D$X) 
# abline(lm(Y9~D$X),col="red")

shade(mu1,D$X,col=col.alpha("blue",0.2))
  

## (2) Linear regression:
# - - - - - - - - - - - -
V = D
colnames(V)  = c("O","B")

V$B = decostand(V$B,"standardize") 

head(V) ; dim(V)  
hist(V$B)
V$O=V$O-2000
plot(V$B~V$O)

bm1 <- quap(
  alist(
    B ~ dnorm(mu, sigma),
    mu ~ a + b*O,
    a   ~ dnorm(0,1),b~dnorm(0,1),sigma ~ dunif (0,50)
  ), data = V )
prec <- precis(bm1,prob=0.95) ; prec

abline(lm(V$B~V$O))

qq = quantile(V$O,probs=c(0,1)) 

weight.seq = seq( from=qq[1] , to=qq[2] , length.out=19 ) 
pred_dat = list( O=weight.seq , O2=weight.seq^2 )

mu = link( bm1 , data=pred_dat ) ; mu.mean <- apply( mu , 2 , mean )
mu.PI = apply( mu , 2 , PI , prob=0.89 ) ; sim.height <- sim( bm1 , data=pred_dat )
height.PI = apply( sim.height , 2 , PI , prob=0.90 )

plot(V$B ~ V$O, col=col.alpha("blue",0.2), 
     xlab="Year", 
     ylab="Total Seed Mass", 
     pch = 16, 
     main="Total Seed Mass",
     cex=1, cex.axis=0.0001, col.axis = "white",
     xlim=c(qq[1],qq[2]),ylim=c(min(V$B)*1.2,max(V$B)*1.2))

lines( weight.seq , mu.mean ,col="blue") 
shade( mu.PI , weight.seq ,col=col.alpha("blue",0.1))
shade( height.PI , weight.seq ,col=col.alpha("blue",0.05))

axis(side=1, at=c(0:23), col="black",labels=c(2000:2023)) 
axis(side=2, at=c(-4,-3,-2,-1,0,1,2,3,4), 
     col="black",
     labels=c(-4,-3,-2,-1,0,1,2,3,4)) 

range(V$B)


#----------------------------------------------------------------------------
## Render fruit and flower production data in phase for the regression models
#----------------------------------------------------------------------------

## Calculate monthly fruit production values:
# - - - - - - - - - - - - - - - - - - - - - -
head(SRM9);dim(SRM9)
SRM9$SRM=SRM2
head(SRM9);dim(SRM9)
#plot(SRM9$SRM,type="l")

SRM8 = SRM9 %>% group_by(MON) %>% summarize(SRM = mean(SRM))
SRM8 = as.data.frame(SRM8)
head(SRM8) ; dim(SRM8)

## Calculate monthly flower production values:
# - - - - - - - - - - - - - - - - - - - - - -
FP9 = SRM9 ; head(FP9) ; dim(FP9)
colnames(FP9)[4] = "TFP"
FP9$TFP = FP2
head(FP9);dim(FP9)

FP8 = FP9 %>% group_by(MON) %>% summarize(TFP = mean(TFP))
FP8 = as.data.frame(FP8)
head(FP8) ; dim(FP8)

## Plot mean monthly values of flower and seed production to decide translation values below:
#--------------------------------------------------------------------------------------------
plot  (SRM8$SRM~SRM8$MON,type="b",col="blue",ylim=c(-30,40))
points(FP8$TFP~FP8$MON,type="b",col="red")
abline(h=0,lty=2)

## Temporal translation of data to render fruit and flower production in phase:
#------------------------------------------------------------------------------
nryrs = length(table(SRM9$YR)[table(SRM9$YR)>0]) ; nryrs
years = names (table(SRM9$YR)[table(SRM9$YR)>0]) ; years

SRM10 = SRM9 ; head(SRM9) ; dim(SRM9)
TFP10 = FP9  ; head(TFP10) ; dim(TFP10)
TFP00 = FP9 
SRM9$YR-FP9$YR ## => ok, years are matching

i=2
for(i in 1:nryrs){
  w                  = which(SRM9$YR%in%years[i])
  SRM10$SRM[w[1:12]] = SRM9$SRM[w[c(9:12,1:8)]]
  TFP10$TFP[w[1:12]] = FP9$TFP [w[c(4:12,1:3)]] ## flower peak 5 months before fruit peak => 5-months translation
  TFP00$TFP[w[1:12]] = FP9$TFP [w[c(9:12,1:8)]] 
}


head(SRM10) ; dim(SRM10) ## Seed Rain Mass data (month-level) in phase with flower data
head(TFP10) ; dim(TFP10) ## Flower data (month-level) in phase with Seed Rain Mass data

FP88 = TFP10 %>% group_by(MON) %>% summarize(TFP = mean(TFP))
FP88 = as.data.frame(FP88)
head(FP88) ; dim(FP88)

FP00 = TFP00 %>% group_by(MON) %>% summarize(TFP = mean(TFP))
FP00 = as.data.frame(FP00)
head(FP00) ; dim(FP00)

SRM88 = SRM10 %>% group_by(MON) %>% summarize(SRM = mean(SRM))
SRM88 = as.data.frame(SRM88)
head(SRM88) ; dim(SRM88)

## Verification that fruit and flower data are in phase:
plot  (SRM88$SRM~SRM88$MON,type="b",col="blue",ylim=c(-30,40))
points(FP88$TFP~FP88$MON,type="b",col="red")
points(FP00$TFP~FP00$MON,type="b",col="orange2")
abline(h=0,lty=2)

head(SRM10);dim(SRM10) ## Community-level seed production data
head(TFP10);dim(TFP10) ## Community-level flower production data

## For each species:
SP.FP1 = log(SP.FP+1) ## FP1 = Flower production 1
head(SP.FP1[,1:10]) ; dim(SP.FP1)

SP.SP1 = log(SP.SRM+1) ## SP1 = Seed production 1
head(SP.SP1[,1:10]) ; dim(SP.SP1)

SP.FP2 = SP.FP1 ; head(SP.FP2) ; dim(SP.FP2)
SP.SP2 = SP.SP1 ; head(SP.SP2) ; dim(SP.SP2)

## We do the same as we did at the community level, i.e. translating
## plots (censuses) for seed production to center the fruit phenology
## within phenological years, and translate the flower production so that the
## latter overlaps with fruits (5 months between peak flowering and fruiting)
nryrs = length(table(SRM10$YR)[table(SRM10$YR)>0]) ; nryrs
years = names (table(SRM10$YR)[table(SRM10$YR)>0]) ; years

i=2
for(i in 1:nryrs){
  w                = which(SRM9$YR%in%years[i])
  SP.FP2[w[1:12],] = SP.FP1[w[c(4:12,1:3)],]
  SP.SP2[w[1:12],] = SP.SP1[w[c(9:12,1:8)],] 
}

head(SP.FP2) ; dim(SP.FP2)
head(SP.SP2) ; dim(SP.SP2)

## verification that translations worked:
SP.FP1[c(1:12),1]
SP.FP2[c(1:12),1] # => ok


######################################################
#*****************************************************
###               Regression Models               ###
#*****************************************************
######################################################

## Seed production data at the year level
#----------------------------------------
## SRM10 and TFP10 are the seed and flower production (community-level)
## in phase at the month level 
## SP.SP2 is the seed production in phase but for each species
dim(SRM10);dim(SP.SP2)
  
## Total fruit production:
SRM10$YR
years = as.numeric(names(table(SRM10$YR))) ; years
nyrs  = length(years) ; nyrs

## Create matrix FR.YR = Seed production at the year level for each species
FR.YR = as.data.frame(matrix(ncol=ncol(SP.SP2),nrow=nyrs))
colnames(FR.YR) = colnames(SP.SP2)
rownames(FR.YR) = years
head(FR.YR) ; dim(FR.YR)

i=2
for(i in 1:nyrs){
  w = which(SRM10$YR==years[i]) ; w
  FR.YR[i,] = as.numeric(colSums(SP.SP2[w,]))
}
head(FR.YR) ; dim(FR.YR)

FR.YR2 = decostand(FR.YR,"standardize")
head(FR.YR2) ; dim(FR.YR2)

FR9 = as.numeric(rowSums(FR.YR2)) ; FR9

## Flower production data at the year level
#------------------------------------------
head(SP.FP2) ; dim(SP.FP2) ; dim(SRM10)
SRM10$YR
years = as.numeric(names(table(SRM10$YR))) ; years
nyrs  = length(years) ; nyrs

FL.YR = as.data.frame(matrix(ncol=ncol(SP.FP2),nrow=nyrs))
colnames(FL.YR) = colnames(SP.FP2)
rownames(FL.YR) = years
head(FL.YR) ; dim(FL.YR)

i=2
for(i in 1:nyrs){
  w = which(SRM10$YR==years[i]) ; w
  FL.YR[i,] = as.numeric(colSums(SP.FP2[w,]))
}
head(FL.YR) ; dim(FL.YR)

FL.YR2 = decostand(FL.YR,"standardize")
head(FL.YR2) ; dim(FL.YR2)

FL9 = as.numeric(rowSums(FL.YR2)) ; FL9
plot(FL9,type="b")

FL9.FR9 = as.data.frame(cbind(FL9,FR9))
colnames(FL9.FR9) = c("FL9","FR9") ; rownames(FL9.FR9) = years
head(FL9.FR9) ; dim(FL9.FR9) 

nryrs = length(table(SRM10$YR)) ; nryrs
years = names(table(SRM10$YR))  ; years


#---------------
## Climate data
#---------------
clim2 <- as.data.frame(read.table("SuppInfo_Appendix_S5.xls",h=T))
head(clim2);dim(clim2)

clim2 = clim2[clim2$YR%in%years,] ; head(clim2) ; dim(clim2)
for(i in 1:nryrs){
  w               = which(clim2$YR%in%years[i])
  clim2[w[1:12],] = clim2[w[c(9:12,1:8)],] 
}

head(clim2);dim(clim2)
clim2$YR

clim3 = clim2[,which(colnames(clim2)%in%c("IRRA","RAIN","TMIN","TMAX","RHAVE"))]
head(clim3);dim(clim3)

E0 = clim3 ## store un-transformed climate data
clim4 = clim3

## Normalize variables:
for(i in 1:ncol(clim4)){
  lambda1 <- BoxCox.lambda(clim4[,i]+abs(min(clim4[,i]))+0.0000000001,method = "loglik", lower = -3,  upper = 3) 
  clim4[,i] <- BoxCox(clim4[,i], lambda=lambda1)
}
par(mfrow=c(3,3),mex=0.3)
for(j in 1:ncol(clim4)){hist(clim4[,j],col="blue",main = colnames(clim4)[j])}
par(mfrow=c(1,1))

## Standardize variables:
for(j in 1:ncol(clim4)){ clim4[,j] <- decostand(clim4[,j],"standardize") }
par(mfrow=c(4,3),mex=0.3)
for(j in 1:ncol(clim4)){hist(clim4[,j],col="lightblue",main = colnames(clim4)[j])}
par(mfrow=c(1,1))
head(clim4);dim(clim4)

## To check correlation between climate variables during certain months:
yrs9 = clim2$YR [clim2$YR%in%years] ; yrs9
mon9 = clim2$MON[clim2$YR%in%years] ; mon9
clim5 = cbind(yrs9,mon9,clim4) ; colnames(clim5) = c("YR","MN",colnames(clim4))
head(clim5);dim(clim5)

clim6 = clim5[clim5$MN%in%c(1,2),] ; head(clim6) ; dim(clim6)

clim7 <- clim6 %>% group_by(YR) %>% summarize(irra  = mean(IRRA),
                                              rain  = mean(RAIN),
                                              tmin  = mean(TMIN),
                                              tmax  = mean(TMAX),
                                              rhave = mean(RHAVE))
clim7 = as.data.frame(clim7)
head(clim7) ; dim(clim7)


#-------------
## FIGURE 3 ##
#-------------
head(SRM10) ; dim(SRM10)

par(mfrow=c(3,2),mex=0.9)

for(i in 1:6){

if(i==1) { Y1 = decostand(as.numeric(SRM10$SRM),"standardize"); tit1 = "(a) Seed Production" ; tity="Seed Production"}
if(i==2) { Y1 = as.numeric(clim2$IRRA)  ; tit1 = "(b) Irradiance" ; tity = "Irradiance (W/m²)"}
if(i==3) { Y1 = as.numeric(clim2$RAIN)  ; tit1 = "(c) Rainfall"   ; tity = "Rainfall (mm)"}
if(i==4) { Y1 = as.numeric(clim2$TMIN)  ; tit1 = "(d) TMIN"       ; tity = "TMIN (°C)"}
if(i==5) { Y1 = as.numeric(clim2$TMAX)  ; tit1 = "(e) TMAX"       ; tity = "TMAX (°C)"}
if(i==6) { Y1 = as.numeric(clim2$RHAVE) ; tit1 = "(f) RHAVE"      ; tity = "RHAVE (%)"}

my.vector3 <- ts(1:240, frequency = 12, start = c(2000, 1)) 
my.vector3[1:240] <- Y1
my.vector3
  
stl1 <- stl(my.vector3, "per") # plot(stl1)
season <- as.numeric(stl1$time.series[,1]) # plot(season,type="l")
intera <- as.numeric(stl1$time.series[,2]) # plot(intera,type="l")

plot  (Y1~SRM10$YRMON,type="l",col="grey70",xlab="Year",ylab=tity, main=tit1,
       cex.lab=1.2)
points(intera~SRM10$YRMON,type="l",col="black")

segments(2017.96,Y1[216],2019.08,Y1[217],col="white",lwd=2,lty=2)
segments(2019.96,Y1[228],2021.08,Y1[229],col="white",lwd=2,lty=2)

#segments(2017.96,Y1[216],2019.08,Y1[217],col="grey60",lwd=1,lty=3)
#segments(2019.96,Y1[228],2021.08,Y1[229],col="grey60",lwd=1,lty=3)

abline(v=seq(from=2018.01,to=2018.99,by=0.01),col="black")
abline(v=seq(from=2020.01,to=2020.99,by=0.01),col="black")

abline(lm(Y1~SRM10$YRMON),lty=2)

}


#-----------------------------------------
## Calculate MEM and MSR for climate data:
#-----------------------------------------
yrmon9 = clim2$YR+(rep(c(1:12),nryrs)/12.5) ; yrmon9
C1 = as.data.frame(cbind(yrmon9,1)) ; head(C1) ; dim(C1) ; plot(yrmon9)
E1 = as.data.frame(clim4) ; head(E1) ; dim(E1)

candidates <- listw.candidates(coord=C1,style="B",nb=c("gab","mst"),weights=c("flin","up"),y_fup = 0.5)
spatial.effect <- 0  
MEM <- listw.select(E1, candidates, MEM.autocor = "positive",
                    method = "FWD", MEM.all = TRUE, nperm = 999,
                    nperm.global = 999, alpha = 0.01, p.adjust = TRUE)      
DN <- 0
if(is.list(MEM)){
  if(is.list(MEM$best)){
    MEM.find           <- 1
    MEM.select         <- as.matrix(MEM$best$MEM.select)
    spatial.effect     <- 1
    DN                 <- candidates[[MEM$best.id]]
    if(i==1){DN2 <- DN}
  }
}
w1 <- which(MEM$best$summary[,6]<=0.001)
MEMs <- MEM.select[,w1]
head(MEMs);dim(MEMs)

## MSR randomization:
nMSR=999
MSR <- msr(E1, DN, nrepet = nMSR, method = "singleton")
MSR[[1]]
E0 <- as.data.frame(E0) ; head(E0) ; dim(E0)
E1 <- as.data.frame(E1) ; head(E1) ; dim(E1)

## Example of MSR-randomized climate data (irradiance) compared with observed values:
par(mfrow=c(2,1),mex=0.75)
plot(E1$IRRA     ~yrmon9,type="l",xlab="year",ylab="Observed solar irradiance")
plot(MSR[[5]][,1]~yrmon9,type="l",col="red",xlab="year",ylab="Randomized solar irradiance")
par(mfrow=c(1,1))

#----------------------
## Regression models ##
#----------------------

## Matrices to store results:
#----------------------------
Max.translations = 12  ## Max number of random translation lag samples
Max.range        = 12  ## Max range (in month) within climate values are calculated

nyears0 = length(table(clim2$YR)[table(clim2$YR)>0]) ; nyears0
years0  = names (table(clim2$YR)[table(clim2$YR)>0]) ; years0

nyears  = length(table(SRM10$YR)[table(SRM10$YR)>0]) ; nyears
years   = names(table(SRM10$YR)[table(SRM10$YR)>0]) ; years

## Matrix B and E to store, using rbind: r2, t2 and correlation values (3 columns)
t2.vec = rep(0,13) ; for(i in 1:12){ t2.vec = c(t2.vec,rep(i,13))} ; t2.vec
r2.vec = rep(c(0:12),13) ; r2.vec
t2r2 = data.frame(t2=t2.vec,r2=r2.vec)
head(t2r2);dim(t2r2)
NR = nrow(t2r2) ; NR

B = as.data.frame(matrix(0,ncol=16,nrow=NR))
colnames(B) = c("r2","t2","A","I","R","T1","T2","H","F",
                "A.sign","I.sign","R.sign","T1.sign","T2.sign","H.sign","F.sign") 
head(B) ; dim(B)

B2 = as.data.frame(matrix(0,ncol=10,nrow=NR))
colnames(B2) = c("r2","t2","A","I","R","T1","T2","H","F","BIC")

B.irra = B2 ; B.rain = B2 ; B.tmin = B2
B.tmax = B2 ; B.have = B2 ; B.flow = B2

B.all = as.data.frame(matrix(0,ncol=20,nrow=NR))
colnames(B.all) = c("BIC","R2","I","R","T1","T2","H","F",
                    "I.t2","R.t2","T1.t2","T2.t2","H.t2","F.t2",
                    "I.r2","R.r2","T1.r2","T2.r2","H.r2","F.r2")

## To test direct and indirect effects via mediator (flower):
med.irra = as.data.frame(matrix(ncol=8,nrow=NR))
colnames(med.irra) = c("t2","r2",
                       "ACME","ADE","TOT",
                       "P.ACME","P.ADE","P.TOT")
head(med.irra);dim(med.irra)
med.rain=med.irra
med.tmin=med.irra
med.tmax=med.irra
med.rhave=med.irra

FL9 = as.numeric(FL9.FR9$FL9) ; FL9
FR9 = as.numeric(FL9.FR9$FR9) ; FR9

## Use year-level sums of climate variables, so render everything positive
## at the census level:
X.CLIM3 = clim4
for(z in 1:ncol(X.CLIM3)){ X.CLIM3[,z] = X.CLIM3[,z] + abs(min(X.CLIM3[,z])) }
head(X.CLIM3) ; dim(X.CLIM3)

## MSR test loop
#---------------
par(mfrow=c(1,1))
plot(1:10,col="white",xlim=c(0,10))

nMSR=4999 ## Number of Moran Spectral Randomizations

R2.I=0 ; R2.R=0 ; R2.T1=0 ; R2.T2=0 ; R2.H=0 ; R2.F=0
X.list = list() ## list object to store X data obtained for each value of t2 and r2
plot(0:10,col="white",xlim=c(0,10),xlab="Loop progression")

for(t in 1:NR){ 
  
  t2 = t2r2$t2[t]
  r2 = t2r2$r2[t]
  B$t2[t] = t2 ; B$r2[t] = r2    ## B stores correlation analyses
  if(r2==12) {vec1=c((24-t2):(36-t2))}
  if(r2<12)  {vec1=c((24-t2+(12-r2)+0):(36-t2))}
  if(r2>12)  {vec1=c((24-t2-r2+12+0):(36-t2))}
  vec2=vec1 ; vec3=vec1 ; vec4=vec1 ; vec5=vec1 ; vec6=vec1
  med.irra$t2[t]  = t2 ; med.irra$r2[t] = r2
  med.rain$t2[t]  = t2 ; med.rain$r2[t] = r2
  med.tmin$t2[t]  = t2 ; med.tmin$r2[t] = r2
  med.tmax$t2[t]  = t2 ; med.tmax$r2[t] = r2
  med.rhave$t2[t] = t2 ; med.rhave$r2[t] = r2
  
  ## Empty year-level Matrix M to be filled:
  tit1 = c("FRUIT","FLOW",
           "IRRA","RAIN","TMIN","TMAX","RHAVE") ; tit1
  M           <- as.data.frame(matrix(ncol=length(tit1),nrow=nyears-1))
  colnames(M) <- tit1
  rownames(M) <- years[2:length(years)]
  head(M) ; dim(M)
  
  M$FRUIT = FR9[-which(years==2000)]
  M$FLOW  = FL9[-which(years==2000)]

  ## Matrix CLIM100 (census-level) of climate variables:
  CLIM100 <- X.CLIM3
  head(CLIM100) ; head(clim2) ; dim(CLIM100) ; dim(clim2) ; dim(SRM10) ; dim(M)

  for(g in 2:nyears){
    
    if(g>3){
      CLIM102 = CLIM100[clim2$YR%in%c(years[g-2],years[g-1],years[g]),]
    }
    
    if(years[g]==2001){ CLIM102 = rbind(CLIM100[c(1:12,1:12,13:24),]) }

    ## Explanatory variables (modified => translation t2 and range r2)
    M$IRRA [g-1] = sum(CLIM102$IRRA[vec1])
    M$RAIN [g-1] = sum(CLIM102$RAIN[vec2])
    M$TMIN [g-1] = sum(CLIM102$TMIN[vec3])
    M$TMAX [g-1] = sum(CLIM102$TMAX[vec4])
    M$RHAVE[g-1] = sum(CLIM102$RHAVE[vec5])

  }
  
  for(j in 1:7){  ## for each explanatory variable...
    
    if(j==1) { X10 <- decostand(M[,which(colnames(M)%in%c("IRRA","RAIN","TMIN","TMAX","RHAVE"))],"standardize") }
    if(j==2) { X10 <- decostand(M$IRRA,"standardize") }
    if(j==3) { X10 <- decostand(M$RAIN,"standardize") }
    if(j==4) { X10 <- decostand(M$TMIN,"standardize") }
    if(j==5) { X10 <- decostand(M$TMAX,"standardize") }
    if(j==6) { X10 <- decostand(M$RHAVE,"standardize") }
    if(j==7) { X10 <- decostand(M$FLOW,"standardize") }
    
    Y10 = decostand(M$FRUIT,"standardize")
    Y10 = as.numeric(decostand(Y10,"standardize"))

    # - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - -
    ## Now calculate the OBSERVED correlation between fruit production and climate:
    # - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - -

    if(j==1){ 
      rda0 = rda(Y10,X10)  
      K = RsquareAdj(rda0)$adj.r.squared 
    }
    if(j>1) { 
        K = lm(Y10~X10)$coefficients[2] 
        R2 = K^2
        if(j==2){if(R2>R2.I) {Y.I  = Y10 ; X.I  = X10 ; R2.I  = R2 ; t.I = t2 ; r.I = r2}}
        if(j==3){if(R2>R2.R) {Y.R  = Y10 ; X.R  = X10 ; R2.R  = R2 ; t.R = t2 ; r.R = r2}}
        if(j==4){if(R2>R2.T1){Y.T1 = Y10 ; X.T1 = X10 ; R2.T1 = R2 ; t.T1= t2 ; r.T1= r2}}
        if(j==5){if(R2>R2.T2){Y.T2 = Y10 ; X.T2 = X10 ; R2.T2 = R2 ; t.T2= t2 ; r.T2= r2}}
        if(j==6){if(R2>R2.H) {Y.H  = Y10 ; X.H  = X10 ; R2.H  = R2 ; t.H = t2 ; r.H = r2}}
        if(j==7){if(R2>R2.F) {Y.F  = Y10 ; X.F  = X10 ; R2.F  = R2 ; t.F = t2 ; r.F = r2}}
      
        if(j<7){
          df            = data.frame(Y=Y10,X=X10,F=decostand(M$FLOW,"standardize"))
          mediate_model = lm(F~X,data=df) ; mediate_model
          full_model    = lm(Y~F+X,data=df) ; full_model
          c_path        = lm(Y~X,data=df) ; c_path
          med1 = mediate(mediate_model,full_model,sims=1,boot=TRUE,treat="X",mediator="F")
          ACME = med1$d0 ; TOT  = c_path$coefficients[2] ; ADE  = TOT-ACME 
          if(j==2){ med.irra$ACME[t]=ACME ;med.irra$TOT[t]=TOT ;med.irra$ADE[t]=ADE}
          if(j==3){ med.rain$ACME[t]=ACME ;med.rain$TOT[t]=TOT ;med.rain$ADE[t]=ADE}
          if(j==4){ med.tmin$ACME[t]=ACME ;med.tmin$TOT[t]=TOT ;med.tmin$ADE[t]=ADE}
          if(j==5){ med.tmax$ACME[t]=ACME ;med.tmax$TOT[t]=TOT ;med.tmax$ADE[t]=ADE}
          if(j==6){ med.rhave$ACME[t]=ACME;med.rhave$TOT[t]=TOT;med.rhave$ADE[t]=ADE}
        }
        
      }
      
      if(j==1) { B$A[t]  = K }
      if(j==2) { B$I[t]  = K }
      if(j==3) { B$R[t]  = K }
      if(j==4) { B$T1[t] = K }
      if(j==5) { B$T2[t] = K }
      if(j==6) { B$H[t]  = K }
    
    ## To store M:
    X.list[[t]] = as.data.frame(M[,which(colnames(M)%in%c("IRRA","RAIN","TMIN","TMAX","RHAVE"))])
    
    null.stat=c() ## vector to be filled with null slope value for the total effect
    null.tot=c()  ## idem (total effect)
    null.acme=c() ## vector to be filled with null slope value for the mediated effect
    null.ade=c()  ## vector to be filled with null slope value for the direct effect

    for(k in 1:nMSR){ 
          
          X.CLIM.null = MSR[[k]] ; head(X.CLIM.null) ; dim(X.CLIM.null)
          X.CLIM.null = as.data.frame(X.CLIM.null)

          M2 = M 
          
          for(z in 1:ncol(X.CLIM.null)){ X.CLIM.null[,z] <- X.CLIM.null[,z] + abs(min(X.CLIM.null[,z])) }

          for(g in 2:nyears){
            
            if(g>3){
              X.CLIM.null2 = X.CLIM.null[clim2$YR%in%c(years[g-2],years[g-1],years[g]),]
            }
            
            if(years[g]==2001){ X.CLIM.null2 = rbind(X.CLIM.null[c(1:12,1:12,13:24),]) }
            
            X.CLIM.null2 = as.data.frame(X.CLIM.null2)
            M2$IRRA [g-1] = sum(X.CLIM.null2$IRRA[vec1])
            M2$RAIN [g-1] = sum(X.CLIM.null2$RAIN[vec2])
            M2$TMIN [g-1] = sum(X.CLIM.null2$TMIN[vec3])
            M2$TMAX [g-1] = sum(X.CLIM.null2$TMAX[vec4])
            M2$RHAVE[g-1] = sum(X.CLIM.null2$RHAVE[vec5])
            
          }
          
          if(j==1) { X1000 = decostand(M2[,which(colnames(M)%in%c("IRRA","RAIN","TMIN","TMAX","RHAVE"))],"standardize") }
          if(j==2) { X1000 = decostand(M2$IRRA,"standardize") }
          if(j==3) { X1000 = decostand(M2$RAIN,"standardize") }
          if(j==4) { X1000 = decostand(M2$TMIN,"standardize") }
          if(j==5) { X1000 = decostand(M2$TMAX,"standardize") }
          if(j==6) { X1000 = decostand(M2$RHAVE,"standardize") }
          if(j==7) { X1000 = decostand(M2$FLOW,"standardize") }
          
          ## Now calculate the OBSERVED correlation between fruit production and climate:
          # - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - -
          if(j==1){ 
            rda1 = rda(Y10,X1000)  
            L = RsquareAdj(rda1)$adj.r.squared 
            null.stat[k] = L
          }
          
          if(j>1) { 
            L = lm(Y10~X1000)$coefficients[2] 
            null.stat[k] = L
            
            if(j<7){
              df            = data.frame(Y=Y10,X=X1000,F=decostand(M$FLOW,"standardize"))
              mediate_model = lm(F~X,data=df) ; mediate_model
              full_model    = lm(Y~F+X,data=df) ; full_model
              c_path        = lm(Y~X,data=df) ; c_path
              med1 = mediate(mediate_model,full_model,sims=1,boot=TRUE,treat="X",mediator="F")
              ACME2 = med1$d0 ; TOT2  = c_path$coefficients[2] ; ADE2  = TOT2-ACME2 
              null.acme[k] = ACME2 ; null.tot[k] = TOT2 ; null.ade[k] = ADE2
            }
            
          }
         
    } ## END "for(k in 1:nMSR){"
        
    ## P-values:
    vec.p = as.numeric(c(K,null.stat)) ; vec.p 
    if(K>0){sig=length(vec.p[vec.p>=as.numeric(K)])/(nMSR+1)}
    if(K<0){sig=length(vec.p[vec.p<=as.numeric(K)])/(nMSR+1)}
        
    if(j==1) { B$A.sign[t]  = sig }
    if(j==2) { B$I.sign[t]  = sig }
    if(j==3) { B$R.sign[t]  = sig }
    if(j==4) { B$T1.sign[t] = sig }
    if(j==5) { B$T2.sign[t] = sig }
    if(j==6) { B$H.sign[t]  = sig }
    if(j==7) { B$F.sign[t]  = sig }
        
    if(j>1){
        
    ## P-values for the direct and indirect effects:
    vec.acme = as.numeric(c(ACME,null.acme)) # vec.acme 
    if(ACME>0){sig1=length(vec.acme[vec.acme>=as.numeric(ACME)])/(nMSR+1)}
    if(ACME<0){sig1=length(vec.acme[vec.acme<=as.numeric(ACME)])/(nMSR+1)}
        
    vec.tot = as.numeric(c(TOT,null.tot)) # vec.tot 
    if(TOT>0){sig2=length(vec.tot[vec.tot>=as.numeric(TOT)])/(nMSR+1)}
    if(TOT<0){sig2=length(vec.tot[vec.tot<=as.numeric(TOT)])/(nMSR+1)}
        
    vec.ade = as.numeric(c(ADE,null.ade)) # vec.ade 
    if(ADE>0){sig3=length(vec.ade[vec.ade>=as.numeric(ADE)])/(nMSR+1)}
    if(ADE<0){sig3=length(vec.ade[vec.ade<=as.numeric(ADE)])/(nMSR+1)}
        
    if(j==2){med.irra$P.ACME[t] =sig1;med.irra$P.TOT[t] =sig2 ; med.irra$P.ADE[t] =sig3}
    if(j==3){med.rain$P.ACME[t] =sig1;med.rain$P.TOT[t] =sig2 ; med.rain$P.ADE[t] =sig3}
    if(j==4){med.tmin$P.ACME[t] =sig1;med.tmin$P.TOT[t] =sig2 ; med.tmin$P.ADE[t] =sig3}
    if(j==5){med.tmax$P.ACME[t] =sig1;med.tmax$P.TOT[t] =sig2 ; med.tmax$P.ADE[t] =sig3}
    if(j==6){med.rhave$P.ACME[t]=sig1;med.rhave$P.TOT[t]=sig2 ; med.rhave$P.ADE[t]=sig3}
        
    }
        
  } ## END "for(j in 1:7)"
  
  points(10/NR*t,5,pch=15,col=rainbow(NR)[t],cex=1)
  
}   ## END "for(t in 1:NR){"


head(B);dim(B) ## total effect of each variable
head(med.irra) ## total, direct and mediated effect of each variable

## climate values from best models for each variable (from next section just below):
XX = as.data.frame(read.table("XX.txt",h=T)) ; head(XX);dim(XX)
YY = as.data.frame(read.table("YY.txt",h=T)) ; YY = as.numeric(YY[,1]) ; YY

med.tmin [which.min(med.tmin$TOT),]
med.rhave[which.min(med.rhave$ADE),]
med.rain [which.min(med.rain$TOT),]
med.rain [which.max(med.rain$ADE),]
med.tmax [which.max(med.tmax$TOT),]
med.irra [which.max(med.irra$TOT),]

#---------------------------------------
###           VIEW RESULTS           ###
#----------------------------------------

### Heatmaps of slope values for each t2 and r2 value:
#-----------------------------------------------------
response  
head(B);dim(B) ## Nr traps and CWM results

#B = as.data.frame(read.table("B.SRM.txt",h=T))
head(B);dim(B) 

H=B ## Y = Nr traps
#H=E ## Y = trait data (FCA)

y=H$r2 ; y
x=H$t2 ; x
n5 = 96
max1 = max(abs(H[,4:8]),na.rm=TRUE) ; max1
min(H$A);H$A[H$A<0]=0.0001;max(H$A)
n6 = max1*1.05 ; n6
head(H);dim(H) 
test.MSR=0

Q <- as.data.frame(matrix(ncol=3,nrow=nrow(H)))
colnames(Q) <- c("x","y","z")
Q$x=x ; Q$y=y 

test.MSR = 1
par(mfrow=c(3,3),mex=0.5)

for(i in 4:8){
  Q$z=as.numeric(H[,i])
  if(i==3){title.graph="All climate variables" 
  if(test.MSR==1){wsig05 = which(H$A.sign <=0.05)}
  if(test.MSR==1){wsig01 = which(H$A.sign <=0.01)}
  if(test.MSR==1){wsig00 = which(H$A.sign <=0.001)}
  }
  if(i==4){title.graph="Irradiance"    
  if(test.MSR==1){wsig05 = which(H$I.sign <=0.05) }
  if(test.MSR==1){wsig01 = which(H$I.sign <=0.01) }
  if(test.MSR==1){wsig00 = which(H$I.sign <=0.001) }
  }
  if(i==5){title.graph="Rainfall"      
  if(test.MSR==1){wsig05 = which(H$R.sign <=0.05) }
  if(test.MSR==1){wsig01 = which(H$R.sign <=0.015) }
  if(test.MSR==1){wsig00 = which(H$R.sign <=0.001) }
  }
  if(i==6){title.graph="TMIN"          
  if(test.MSR==1){wsig05 = which(H$T1.sign<=0.05) }
  if(test.MSR==1){wsig01 = which(H$T1.sign<=0.01) }
  if(test.MSR==1){wsig00 = which(H$T1.sign<=0.001) }
  }
  if(i==7){title.graph="TMAX"          
  if(test.MSR==1){wsig05 = which(H$T2.sign<=0.05) }
  if(test.MSR==1){wsig01 = which(H$T2.sign<=0.01) }
  if(test.MSR==1){wsig00 = which(H$T2.sign<=0.001) }
  }
  if(i==8){title.graph="RHAVE"          
  if(test.MSR==1){wsig05 = which(H$H.sign <=0.05) }
  if(test.MSR==1){wsig01 = which(H$H.sign <=0.01) }
  if(test.MSR==1){wsig00 = which(H$H.sign <=0.001) }
  }
  if(i==9){title.graph="Flowers"       
  if(test.MSR==1){wsig05 = which(H$F.sign <=0.05) }
  if(test.MSR==1){wsig01 = which(H$F.sign <=0.01) }
  if(test.MSR==1){wsig00 = which(H$F.sign <=0.001) }
  }
  mba.int <- mba.surf(Q, n5, n5, extend=T)$xyz.est
  if(i==3){fields::image.plot(mba.int, xlab="", ylab="", axes=F, zlim=c(0,max(H$A*1.05)))}
  if(i>3) {fields::image.plot(mba.int, xlab="", ylab="", axes=F, zlim=c(-n6,n6))}
  #axis(1, at=c(0:12),las=1,cex.axis=1,pos=0)
  #axis(2, at=c(0:12),las=1,cex.axis=1,pos=0)
  title(xlab="Translation lag (in months)", line=0.5, cex.lab=1.2)
  title(ylab="Range (in month)", line=0.5, cex.lab=1.2)
  title(main=title.graph, line=0.8, cex=2)
  #if(test.MSR==1){points(Q$x[wsig05],Q$y[wsig05],pch=15,col="green",cex=1.0)}
  if(test.MSR==1){
    if(length(wsig01)>0){points(Q$x[wsig01],Q$y[wsig01],pch=4,
                                col=col.alpha("black",1),cex=1)}
  }
  #if(test.MSR==1){points(Q$x[wsig00],Q$y[wsig00],pch=17,col="red",cex=1.0)}
}

## Crosses in the heatmaps indicate p-values < 0.05


#-----------------
##   Figure 4   ##
#-----------------

colfunc1 <- colorRampPalette(c("red","grey95","blue"))
colvec <- colfunc1(100)
plot(1:100,col=colvec)
seq1 = seq(from=-0.99,to=0.99,by=0.02) ; seq1

plot(1,1,cex=10,pch=15,col=colvec[60])
plot(1,1,cex=10,pch=15,col=colvec[40])
plot(1,1,cex=10,pch=15,col=colvec[50])

main3 = "Best Model" ## title for each graph of the last graph column

store.BPR.BNR = as.data.frame(matrix(ncol=6,nrow=5))
colnames(store.BPR.BNR) = c("BPR.TOT","BNR.TOT","BPR.ADE","BNR.ADE","BPR.ACME","BNR.ACME")
rownames(store.BPR.BNR) = c("IRRA","RAIN","TMIN","TMAX","RHAVE")

alpha1=0.05
par(mfrow=c(5,4),mex=0.6)
#plot(1,1) ; plot(1,1)

u=5;h=1;v=1
for(u in 1:5){ ## for each climate variable...

if(u==1){MED=med.irra ; main1 = "IRRA"  ; X9 = XX$IRRA   ; tit1 = "IRRADIANCE"}
if(u==2){MED=med.rain ; main1 = "RAIN"  ; X9 = XX$RAIN   ; tit1 = "RAINFALL"}
if(u==3){MED=med.tmin ; main1 = "TMIN"  ; X9 = XX$TMIN   ; tit1 = "TMIN"}
if(u==4){MED=med.tmax ; main1 = "TMAX"  ; X9 = XX$TMAX   ; tit1 = "TMAX"}
if(u==5){MED=med.rhave; main1 = "RHAVE" ; X9 = XX$RHAVE  ; tit1 = "RHAVE"}

  
for(h in 1:3){ ## for each of the 3 statistics ("TOT","ADE","ACME")

if(h==1){
  diff1 = (24-MED$t2)+(MED$r2/2) ; diff1
  MED2  = MED[which(diff1<=24),] ; head(MED2);dim(MED2) 
  diff2 = (24-MED2$t2)-(MED2$r2/2) ; diff2
  MED2  = MED2[which(diff2>=1),] ; head(MED2);dim(MED2)  
  MED2  = MED2[MED2$r2>=1,]
  MED2  = MED2[-which(MED2$t2<=5 & MED2$r2<=5),]
} 

if(h==2){ 
  diff1 = (24-MED$t2)+(MED$r2/2) ; diff1
  MED2  = MED[which(diff1<=24),] ; head(MED2);dim(MED2) 
  diff2 = (24-MED2$t2)-(MED2$r2/2) ; diff2
  MED2  = MED2[which(diff2>=1),] ; head(MED2);dim(MED2) 
  MED2  = MED2[MED2$r2>=1,]
  MED2  = MED2[-which(MED2$t2<=5 & MED2$r2<=5),]
}
  
if(h==3){
  diff1 = (24-MED$t2)+(MED$r2/2) ; diff1
  MED2  = MED[which(diff1<=19),] ; head(MED2);dim(MED2)
  diff2 = (24-MED2$t2)-(MED2$r2/2) ; diff2
  MED2  = MED2[which(diff2>=1),] ; head(MED2);dim(MED2) 
  if(u==4)           { MED2  = MED2[MED2$r2>=3,] }
  if(u%in%c(1,2,3,5)){ MED2  = MED2[MED2$r2>=2,] }
}  
  
if(h==1){ stat = MED2$TOT  ; main2 = "TOTAL"    ; pvals = MED2$P.TOT  ; u1=1 ; u2=2 }  
if(h==2){ stat = MED2$ADE  ; main2 = "DIRECT"   ; pvals = MED2$P.ADE  ; u1=3 ; u2=4 }
if(h==3){ stat = MED2$ACME ; main2 = "MEDIATED" ; pvals = MED2$P.ACME ; u1=5 ; u2=6 }
  
store.BPR.BNR[u,u1] = stat[which.max(stat)]
store.BPR.BNR[u,u2] = stat[which.min(stat)]

## mean monthly fruit data at community-level
plot(c(13:24),decostand(SRM88$SRM,"standardize"), col="black",type="l",
      lty=1,xlim=c(6,27),ylim=c(-2,2.3),ylab="",xlab="Month",
      col.axis="white",cex=0.6,pch=16,cex.axis=0.8, main=main2)
  
for(v in 1:nrow(MED2)){
  mid = 24-MED2$t2[v] ; mid
  ran = MED2$r2[v]/2  ; ran 
  val = stat[v]      ; val
  col10="grey80"
  if(pvals[v]<=alpha1){col10=colvec[which.min(abs(seq1-val))]}
  if(pvals[v]> alpha1){col10="grey80"}
  segments(x0=mid-ran,y0=val*2,x1=mid+ran,y1=val*2,
          lwd=0.1,col=col10)
}

wm1 = which.max(abs(stat)) ; wm1
mid = 24-MED2$t2[wm1] ; mid
ran = MED2$r2[wm1]/2  ; ran 
val = stat[wm1] ; val
if(pvals[wm1]<=alpha1){
  if(val<0){col10="red"}
  if(val>0){col10="blue"}
  #col10="green"
  segments(x0=mid-ran,y0=val*2,x1=mid+ran,y1=val*2,lwd=2,col=col10)
  abline(v=c(mid-ran,mid+ran),lty=1,col=col10)
}  

col11="grey80"
wm9 = which.min(stat) ; wm9
mid = 24-MED2$t2[wm9] ; mid
ran = MED2$r2[wm9]/2  ; ran 
val9 = stat[wm9] ; val9
verif9 = val*val9 ; verif9
if(pvals[wm9]<=alpha1 & verif9<0){
  if(val9<0){col11="red"}
  if(val9>0){col11="blue"}
  segments(x0=mid-ran,y0=val9*2,x1=mid+ran,y1=val9*2,lwd=2,col=col11)
}

## Mean monthly fruit data at community-level
points(c(13:24),decostand(SRM88$SRM,"standardize"), col="black",type="l",
       lty=1,xlim=c(5,25),ylim=c(-2,2.3),ylab="",xlab="Time (months)",
       col.axis="white",cex=0.6,pch=16,cex.axis=0.8) 
points(c(24,25),decostand(SRM88$SRM,"standardize")[c(12,1)],
       col="black" ,type="l",lty=1,cex=0.6,pch=16,cex.axis=0.8) 
abline(h=0,lty=3,col="black")
axis(side=2,at=c(-2,0,2),labels=c(-1,0,1),cex.axis=0.8)    
axis(side=1,at=c(6:27),labels=c(c(6:27)),cex.axis=0.8,col.axis="white")    

monthnames = c("F","M","A","N","F","M","A","N")
axis(side=1,at=c(6,9,12,15,18,21,24,27),labels=monthnames,cex.axis=0.8)    

## mean monthly flower data at community-level
points(c(13:20),decostand(FP00$TFP,"standardize")[c(1:8)],
       col="grey10",type="l",lty=2,ylab="",xlab="time (months)",
       main=tit1,cex=0.6,pch=16,cex.axis=0.8) 
points(c(9:13),decostand(FP00$TFP,"standardize")[c(9:12,1)],
       col="grey10",type="l",lty=2,ylab="",xlab="time (months)",
       main=tit1,cex=0.6,pch=16,cex.axis=0.8) 
points(c(8:9),decostand(FP00$TFP,"standardize")[c(8,9)],
       col="grey10",type="l",lty=2,ylab="",xlab="time (months)",
       main=tit1,cex=0.6,pch=16,cex.axis=0.8)

} ## END "for(h in 1:3){"
  
w1 = which.max(abs(MED2$TOT)) 

tot1  = MED2$TOT[w1]
acme1 = MED2$ACME[which.max(abs(MED2$ACME))]
ade1  = MED2$ADE [which.max(abs(MED2$ADE))]
ptot1 = MED2$P.TOT[w1]

df9=data.frame(Y9=YY,X9=X9)

plot(Y9~X9,data=df9,col="black",pch=16,cex=0.8,cex.axis=0.8,main=main3,
     xlim=c(-2.5,2.5),ylim=c(-2.5,2.5),ylab="SP",xlab=tit1,col.axis="white")
abline(h=0,lty=3,col="grey50")
axis(side=1,at=c(-2,0,2),labels=c(-1,0,1),cex.axis=0.8)    
axis(side=2,at=c(-2,0,2),labels=c(-1,0,1),cex.axis=0.8)    

bm1 <- quap(
  alist(
    Y9 ~ dnorm(mu, sigma),
    mu ~ a + b*X9,
    a   ~ dnorm(0,1),b~dnorm(0,1),sigma ~ dunif (0,50)
  ), data = df9 )
prec <- precis(bm1,prob=0.95) # prec
weight.seq = seq( from=-2.5 , to=2.5 , length.out=19 ) 
pred_dat = list( X9=weight.seq )

mu = link( bm1 , data=pred_dat ) ; mu.mean <- apply( mu , 2 , mean )
mu.PI = apply( mu , 2 , PI , prob=0.95 ) ; sim.height <- sim( bm1 , data=pred_dat )
height.PI = apply( sim.height , 2 , PI , prob=0.95 )

cor2  = cor(YY,X9) ; cor2
plaus = prec[2,4]*prec[2,3] ; plaus

if(plaus<=0){
  lines( weight.seq , mu.mean ,col="grey30",lty=1)
  shade( mu.PI , weight.seq ,col=colvec[50])
}
if(plaus>0){
  
if(ptot1<=alpha1){
  if(cor2<0){
    lines( weight.seq , mu.mean ,col="red",lty=1)
    #shade( mu.PI , weight.seq ,col=col.alpha("red",0.2))
    shade( mu.PI , weight.seq ,col=colvec[40])
    
  }
  if(cor2>0){
    lines( weight.seq , mu.mean ,col="blue",lty=1)
    #shade( mu.PI , weight.seq ,col=col.alpha("blue",0.2))
    shade( mu.PI , weight.seq ,col=colvec[60])
  }
}
  
if(ptot1>alpha1){
  lines( weight.seq , mu.mean ,col="black",lty=1)
  #shade( mu.PI , weight.seq ,col=col.alpha("lightgrey",0.2))
  shade( mu.PI , weight.seq ,col=colvec[50])
}
  
}

if(u==1){abline(h=0,lty=3,col="black")}

points(Y9~X9,data=df9,col="black",pch=16,cex=0.8,cex.axis=0.8,
       xlim=c(-2.5,2.5),ylim=c(-2.5,2.5),ylab="",xlab="")

for(k in 1:2){  

if(k==1){coef1=acme1}
if(k==2){coef1=ade1}
  
to=2
from=-2
length.out=20
seq5 = seq(from=-2.5,to=2.5,by = ((to - from)/(length.out - 1)))

vecmed=c();for(i in seq5){vecmed=c(vecmed,coef1*i)} ; vecmed
if(coef1<0){col11="red"} ; if(coef1>0){col11="blue"}
#abline(lm(vecmed~c(-2,-1,0,1,2)),col=col11,lwd=2)
if(u==1){col11="black"}

if(k==1){abline(lm(vecmed~seq5),col=col11,lty=2)}
if(k==2){abline(lm(vecmed~seq5),col=col11)}

} ## END "for(k in 1:2){"

if(u==2){
  ade1  = med.rain$ADE[which.min(med.rain$ADE)]
  to = 2
  from = -2
  length.out = 20
  seq5 = seq(from=-2.5,to=2.5,by = ((to - from)/(length.out - 1)))
  vecmed=c();for(i in seq5){vecmed=c(vecmed,ade1*i)} ; vecmed
  abline(lm(vecmed~seq5),col="red")
}

} ## END "for(u in 1:5){"


